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We present a Green's function method for the evaluation of the particle density profile and of the 
higher moments of the one-body density matrix in a mesoscopic system of TV Fermi particles moving 
independently in a linear potential. The usefulness of the method is illustrated by applications to a 
Fermi gas confined in a harmonic potential well, for which we evaluate the momentum flux and kinetic 
energy densities as well as their quantal mean-square fluctuations. We also study some properties of 
the kinetic energy functional Ek. in [n(x)] in the same system. Whereas a local approximation to the 
kinetic energy density yields a multi-valued function, an exact single-valued relationship between 
the density derivative of Eki n [n(x)] and the particle density n(x) is demonstrated and evaluated for 
various values of the number of particles in the system. 
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C3 . I. INTRODUCTION 

The study of ultracold gases of fermionic alkali atoms is being pursued experimentally M by using basically the same 
trapping and cooling techniques which have led in 1995 to the achievement of Bose-Einstein condensation in gases of 
bosonic isotopes . The main focus of the experimental effort is towards the realization of novel superfluids from 
pairing between fermions in different hyperfine states. However, attention has also been given to dilute Fermi gases 
in a fully spin-polarized state under magnetic confinement. Such a system presents special motives of fundamental 
interest from the fact that the Pauli principle suppresses the interactions in the s-wave scattering channel and thereby 
essentially isolates the kinetic energy functional which is invoked in the density functional theory of inhomogeneous 
fluids §. 

I | A quasi-one dimensional Fermi gas can be realized in the experiments [jD inside a very anisotropic magnetic trap, 
^ . producing approximately harmonic confinement which is very tight in two directions. The particle and kinetic energy 
t-H ' densities of the corresponding one-dimensional (ID) model of non-interacting fermions in a harmonic potential well 
have been evaluated exactly up to mesoscopic numbers of particles ||. The shell structure in the particle density 
which was noticed in earlier studies of 3D Fermi gases |7|,|| is greatly enhanced in lower dimensionality (see also 
Ref. H). As a further motive of interest, the ID model of noninteracting fermions has the same particle density 
profile and thermodynamic properties as the ID model of point-like bosons with hard core interactions pO|-p^|. The 
£^ \ system of interacting bosons in ID is attracting attention due to its rich phase diagram Jl3]] . 

The particle and kinetic energy densities are the zero-th moment and a second moment of the one-body density 
matrix. In our earlier study B we have outlined a general method by which one may calculate the successive 
moments of the density matrix in a system of N Fermi particles moving independently in a ID potential. In the fully 
quantum regime alternative definitions can be given for the higher-order moments, depending on how the operators 
are symmetrized [ fl4| . In particular, two different second-order moments can be defined for an inhomogeneous gas 
under confinement, yielding two different physical quantities which are the kinetic energy density and the momentum 
flux density. The method here proposed is suitable for the higher-order moments in any of their forms. It exploits 
the definition of the Green's function operator G(x) in coordinate space and allows one to use the numerical methods 
already developed in solid state physics fl5]-|l7|| to evaluate the density of single-particle states from the Green's 
, function in the energy domain. 

In this paper we give a detailed account of the general method in Sec. || and some new applications to the case of 
harmonic confinement in the following Sections. In Sec. pTJ| we contrast the kinetic energy density with the momentum 
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flux density in the presence of inhomogeneity. In Sec. IV we assess the quantal (zero-temperature) fluctuations in the 
kinetic energy and momentum flux densities by means of a calculation of fourth moments of the one-body density 
matrix. Properties of the kinetic energy functional are discussed in Sec. |y| in the context of density functional theory. 
An exact local-density method exists for the system of present interest, as shown in a different approach by Lawes and 
March Go], and is explicitely exhibited here for various values of the number of Fermi particles. Finally, in Sec 
we summarize our main conclusions. 



VI 



1 



II. THE METHOD 



We consider a system of N non-interacting Fermi particles in one spatial dimension subjected to an external 
confinement given by a real potential V ex t(x). We may take the eigenf unctions ipi( x ) as rea l [l9| - The one-body 
density matrix at zero temperature is 



p(x, Xi) = }^ ipi(x)ipi(xi). 



(1) 



Its zero-th moment is the equilibrium density profile n(x), and its odd moments are zero due to the reality of the 
wavefunctions. For the even moments of higher order several definitions have been considered in the literature (for a 
review see e.g. Ziff et al. [fblf ), among which we shall focus on the following: 



P n (x) = {-if 



and 



S n (x) = (-i) n 



Qn 

Or 



-p(x + r/2,x-r/2) 



(2) 



(3) 



r=0 



(we have set h — 1 and m = 1). For n > 2 these definitions give rise to different equilibrium properties in an 
inhomogeneous system, although their semiclassical limits and their integrals coincide. In the specific case of n — 2 
Eqs. (g) and (|J) give twice the kinetic energy density and the momentum flux density, respectively: these will be 
discussed in Sec. III . 

The method developed in Ref. Q is here extended to evaluate both P n (x) and S n (x). For this purpose we write 
Eqs. (Q) and (||) in terms of the Green's function G(x) — (x — x + ie)^ 1 in coordinate space as 



1 N 

P n {x) = lim Im \G(x)f l \ ft) 

7T e^0+ * — ' 



(4) 



and 



1 n ( \ N 

= - 2^ £ u ) A Im ?>* \p a G{x)p n ~' j \ii). 



<T = 



(5) 



Here, | tpi) are the eigenstates of the Hamiltonian in the coordinate representation and p is the momentum operator. 
Equation (Q) is easily proved by making use of the expression of G(x) in the coordinate representation: 
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P n (x) = lim Im} dxiipi(xi) 



7T e-»0+ 
JV 



f 



d" 



(-i) n y^ j J dxitpi(xi)6(x - Xi)-g-^tpi(xi) , 



(6) 



where Xi are the eigenvalues of the position operator x. Equation (^) can be proved in a similar fashion. 

The momentum distribution n{p) can be analogously written in terms of the trace of the Green's function G(p) = 
(p—p+is)^ 1 in momentum space, taken on the eigenvectors | (pi) of the Hamiltonian in the momentum representation: 



JV 



(7) 



or 



n{p) 



1 N 

- lim Im yUi\G(p)\ 



(8) 
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where pi are the eigenvalues of the momentum operator. 

To evaluate the Green's functions in the specific case of harmonic confinement, we make use of the representation for 
the position and the momentum operators in the basis of the eigenstates of the harmonic oscillator, i.e. x = (a+a)) / 
and p = i(a) — a)/^/2 with a \ ipi) = \Ji — 1 1 ipi-i) an d | ipi) — Vi\ipi+i) (we have set the harmonic oscillator 
frequency w = 1). Explicitly, the representation of x in matrix form is given by 



/ l/y/2 

1/V2 



V 



1 



^3/2 



V2 
y/2 



(9) 



and similarly for p. The various moments of the one-body density matrix and the momentum distribution are then 
evaluated as the trace of the operators p° 'G(x)p n ~ cr and of G(p) on the first N states (Trjy)- The trace is calculated 
by the general method given below. The quantities n(x) and n{p) could also be calculated from a suitably modified 
Kirkman-Pendry relation [J2 0|J 1 7j] , as is shown in Appendix A. 



A. Expression for the evaluation of the trace 

The partial trace Trjv of a generic matrix Q is related to the determinant of the inverse matrix Q^ 1 by 

Tr N Q = d [lndet(Q- 1 + AIjv)] /3A| A=0 . (fO) 
This relation is demonstrated as follows: 



d [lndet(g- 1 + Xl N )] /9A| A=0 = d [Trln^" 1 + \l N )] /9A| A=0 
= Tr [(Q- 1 + AIat) _1 Ijv] a=0 = Tr(QI N ) = Tr^Q. 



(11) 



Here we take Q = p a G(x)p n a . In the case of harmonic confinement this matrix is the product of a (2n + l)-diagonal 
matrix and of the inverse of a tridiagonal matrix. Exploiting Eq. ( |To| ) we can write 



Pn(x) 



-— lim Im— 

7T e^0+ OA 



lndet(x + ie — K 



r^,0^ 



A=0 



and 



S n (x) = — - — | n | lim Im-^- 



\ndet{x + ie - K n ' a ) 



x=o 



(12) 



(13) 



where k n > a = x — p n ~ a XlNP a ■ We have thus reduced the problem to the evaluation of the determinant of the matrix 
(x + ie — K n ' a ), which is tridiagonal in the tail and whose first N rows have only a few non-vanishing elements for 
low values of n. 

If we make a partition in blocks of dimension n x n for the first part of the matrix, the whole operator can be 
written as 



/ Ai Bi, a 

Z?2,l A2 #2,3 



K n ' a = lim K n A f = lim 



B. 



3,2 



B. 



3,4 



Bm-\m 
Bm,m-i Am 



(14) 



Here, we have introduced the truncated matrix K^f , which will be used below in actual numerical calculations. 
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B. The renormalization procedure 



The determinant of the matrix (x + ie — K n ' a ) with a tridiagonal representation such as that in Eq. ( [14|) can be 
factorized into a product of determinants of matrices having the same dimension as that of the blocks of the partition 

As a first step, it is easy to show that if we partition the matrix K M a into four blocks, thereby defining Bij an 
K™' a from 

Kf = ( Bm /' m ) , (15) 

we can write 

det(x + ie - K M °) = det(x + ie — K M a _^) ■ det(x + ie - Am — Bm,m-i(% + ie — KmL\)~ &m-i,m)- (16) 
Applying recursively this procedure and taking the limit M — > oo we obtain 

M 

det (x + ie- K n ' a ) = dct (x + ie- Ax) ■ hm TT det (x + ie- Aj-Bj,j-x(x + ie~ Kjj^Bj-ij) . (17) 

It is now important to notice that, owing to the particular form of the matrices Bj-xj and Bj t j-i, it is not necessary 
to explicitly invert the matrices [x + ie — K^). Rather, we may use a renormalization procedure for the operator 

K^Li to further simplify the expression (^) and to calculate the inverse of matrices having dimension at most equal 
to n x n. Renormalization is based on the idea that if one only needs to describe a given subspace of the whole 
Hilbert space of the system, one can define an effective operator acting in the subspace and taking into account the 
contribution from the rest of the system. Renormalization allows us to write 

oo 

det(x + ie - K n '°) = J| det (a; + ie - Aj ) (18) 
j=i 

where Ai — A\ and 

Aj = Aj + B jtj ^i(x + ie- Aj-x^Bj^xj (19) 

for j > 1. The Bjj+x have been introduced in Eq. In the specific case of n = 2 and a = Eq. (|l8|) yields back 

the expression obtained in Ref . ||] . 

Operatively, we have studied the convergence of the function det(x+ie — K M a ) on increasing M. We have performed 
our calculations up to n = 4 and N = 1000 fermions with M ~ 10 7 and e = 0.001. 



III. KINETIC ENERGY AND MOMENTUM FLUX DENSITIES 



The second moments P2(%) and ^(x) of the one-body density matrix from Eqs. (|^) and (^) have different physical 
meanings and show different behaviours. The function i"2(x) = (ipi\8(x — Xi)p 2 \t/ji) is twice the kinetic energy 
density T(x) and for harmonic confinement has already been evaluated in Ref. Q together with the equilibrium 
density profile. This quantity is of main interest in the context of Density Functional Theory, as will be discussed in 
Sec. ^. On the other hand, Sa(x) = (1/2) (tpi\S(x — Xj)fP + pS(x — Xi)p\ipi) is the momentum flux density II(x) 
which enters the equations of generalized hydrodynamics pjj . In 3D fluids this quantity becomes a tensor, known as 
the kinetic stress tensor. 

The equation of motion of the particle density nix, t) for non-interacting fermions in ID is 



d 2 d 2 d 

_ n(M) „n (M ) + _ 



d 

n(x,t) — V ext (x,t) 



(20) 



At equilibrium in a static potential we obtain from Eq. 
and the density profile, 



j) an exact relation between the momentum flux density 



dx 



-n(x)—V ext (x) 
ax 



const 



(21) 
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The integration constant can be fixed by imposing suitable boundary conditions for x — > ±00. 

Another exact relation can be derived by exploiting Eqs. (fj])-([|) together with the definition of n(x) as the zero-th 
moment of the density matrix. This reads: 

U(x)=2T(x) + ^n(x), (22) 

Equation (^2|) exhibits the relation which exists between the momentum flux density and the kinetic pressure P(x) = 
2T{x) in the inhomogeneuos system. Evidently, in the homogeneous limit the momentum flux density and the kinetic 
pressure coincide - and therefore they also coincide within a local-density description for the confined gas. They 
instead differ within a quantum-mechanical description of the the inhomogeneous system, when one keeps exactly 
into account the role of confinement. This difference is readily illustrated in the fully quantum case of a single fermion 
in an harmonic oscillator potential, where we have 

H(a) = T(x) + V(x) (23) 

with V(x) — n(x)V ex t(x) being the potential energy density associated with the confinement. 

For the case of harmonic confinement we have employed the method described in Sec. || to evaluate the exact 
profiles for both T(x) and II(a;) at various numbers of particles. Figure [j] shows these results for N= 4, 12 and 24 
fermions. While T(x) shows N oscillations and negative tails II(x) is everywhere positive and its shell structure 
is visible only in its first derivative (see inset in Fig. ^). 



IV. QUANTUM FLUCTUATIONS OF THE KINETIC ENERGY AND MOMENTUM FLUX DENSITIES 



From the fourth moments Si(x) and Pa{x) of the density matrix we can evaluate the local fluctuations of the 
momentum flux and kinetic energy densities. For the former the fluctuations are most conveniently introduced 
through the Wigner distribution function, 



fw(p, x)= dr e lpr p(x + r/2, x - r/2) 



(24) 



In this representation the relevant moments are given by n{x) = fw(P: x ) an d by H(x) = J2 P P 2 fw(p,x) 
n{x)(p 2 ). The mean square fluctuation of the momentum flux density is therefore given by 



An(x) 



(p 2 - (p 2 )Yfw(p,x) = S 4 {x) - U\x)/n{x) , 



(25) 



where S±{x) is as defined in Eq. 
rewritten as 



Afl(i 



Equation ( p5| ) in the coordinate representation for the density matrix can be 

n(x) " 



dr 2 



n(x) 



p(x + r/2,x-r/2)\ r=0 . 



(26) 



The fluctuations of the kinetic energy density are obtained in a similar fashion as 



AT(x) 



1 d 2 
~2dx\ 



n(xi) 
T 2 {x) | 
n(x) 



1 d 

2dx 



n{x) Tx 



T(x) 
n(x) 



(27) 



with Pi(x) as defined in Eq. (S). 

For harmonic confinement in the local density approximation (LDA) we have tilda{x) = (2N—x 2 ) 1 ^ 2 /n, Hlda(x) = 
2T LDA (x) = (2N-x 2 ) 3 / 2 /3n and S% DA {x) = P^ DA {x) = {2N - x 2 ) 5 / 2 /5n. It is straightforward to show from these 
expressions that the integrated mean square fluctuations, divided by the square of the corresponding integrated 
quantity, scale as 1/N. The same scaling is found in the exact numerical results at large N. 

We have used the Green's function method described in Sec. [n] to evaluate the exact profiles for both ATl(x) and 
AT(x) in harmonic confinement, as shown in Fig. [2] for TV =4, 12 and 24 fermions. The positions of the N — 1 maxima 
in AT(x) are located in correspondence to the minima of the kinetic energy density in Fig. [I] (a), and again negative 
tails are present at the boundaries. In the case of N = 2 we have tested our numerical results against the analytic 
expression given for AT(x) by low-order Hermite polynomials, the result being shown in Fig. |. Evidently, we have 
obtained in this case a very accurate numerical evaluation of up to the fourth moment of the density matrix. 
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V. PROPERTIES OF THE KINETIC ENERGY FUNCTIONAL 



We have seen that the second moment P2(x) of the one-body density matrix yields the kinetic energy density 
T{x) = P2{x)/2. The total kinetic energy is obtained by integration over the spatial coordinate and in the spirit of 
Density-Functional Theory (DFT) f§ can be viewed as a functional of the particle density n(x): 

E k in[n(x)] = { dxT(x) = fdxf[n(x)] . (28) 



Here we have introduced an operator T which acts on the density profile to yield 

f[n(x)} = T(x) . (29) 

A main point of interest for DFT is to assess the properties of the unknown operator T. 

Knowing exactly the function T(x) the specific case of ID harmonic confinement, we may try to construct a function 
T(n) — T(x(n)) in a local density approach, with x(n) being given by inversion of the density profile. The inversion 
has to be performed locally, since the exact density profile is not globally invertible owing to its shell structure. The 
resulting T(n), as reported in Fig^, is a multi- valued function around the values of the particle density corresponding 
to each local maximum in the profile. While this shows that a local-density approach cannot characterize the operator 
T, further progress can be made by studying the derivative of T(n) with respect to the density, which is the local 
chemical potential jlioc{n) = dT(n)/dn. 

We show here that for the specific case of a ID harmonic confinement the inverse n(/x; oc ) of the function fli oc (n) is 
a single-valued function. This is most simply seen by inverting the Euler equation 

5E Hn [n(x)\ _ 

— j — 1-\ — = fJ- - V ext {x) = nioc(x) (30) 
onlyX) 

in the domain x > and by substituting its value into the density profile n{x) to obtain 

^(Woc) = n(x(pi oc )) . (31) 

Here, x(/j,i oc ) is the inverse of the function /x; oc (a;). Because of the symmetry x «-> —x for the harmonic potential, the 
region x < is not needed to obtain n(fii oc ). 

An analytic calculation for n(/i; oc ) has been reported by Lawes and March [jl8| for N = 1 and 2. The numerical 
method described in Sect. || allows us to easily compute the density profile for numbers of particles up to N = 1000. 
By combining these results for n(x) with the analytic expression for x(/i; oc ), we have obtained the function n(ni oc ) 
which is shown in Fig. [| for various values of N. At variance from the LDA prediction, the exact density profile is 
non- vanishing at negative values of [ii oc and shows a number of oscillations (equal to N/2 for even N or to (N — l)/2 
for odd N), with decreasing amplitude as N increases. 

The property of single- valuedness of n(fii oc ) is a characteristic of ID harmonic confinement and can be understood 
to be a consequence of the local nature of the relation n(fii oc ) vs. ni oc , which was already pointed out by Lawes and 
March Anharmonicity would be sufficient to make it invalid. 



VI. SUMMARY AND CONCLUSIONS 



In summary, we have presented a general method for evaluating the n-th moments of the one-body density matrix 
(defined either through the derivatives with respect to the relative coordinate or through the derivatives with respect 
to the second position variable) for a system of confined non-interacting ID fermions. We have applied this method 
to the case of harmonic confinement to calculate up to fourth moments. 

Regarding second moments, the momentum flux density H(x) has been computed and compared with the kinetic 
energy density T(x) at various numbers of particles. While T(x) shows a prominent shell structure and negative tails, 
H(x) is everywhere positive and has a less marked shell structure, which can be evidenced in its first spatial derivative. 
The relation between these two different physical quantities has been elucidated by giving an exact relation between 
them and the equilibrium density profile of the fermion cloud. The momentum flux density determines the dynamical 
properties of the cloud, whereas the kinetic energy density is relevant in the context of Density Functional Theory. 

We have also calculated the fourth moments of the density matrix in order to display the quantal fluctuations of the 
kinetic energy and momentum flux densities. Finally, we have evaluated the local relationship which exists between 
density profile and local chemical potential for fermion clouds in ID harmonic confinement. This relation has allowed 
us to display an exact property of the DFT kinetic energy functional for a ID Fermi gas up to mesoscopic values of 
the number of particles. 
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APPENDIX A: DENSITY PROFILE FROM THE KIRKMAN-PENDRY RELATION 



We derive in this Appendix an expression alternative to Eq. (|12j) for the evaluation of the density profile in a 
system of N fermions under ID harmonic confinement. The trace of the Green's function G(x) over its first N states 
in Eq. (Q) is obtained through an extension of the Kirkman-Pendry relation ]2(J , already used for the density of states 



in a quasi- ID system |15 



1( 



We define the Green's function 

G(S,x) = -4— . (Al) 

X + — £(x) + IE 

where 5 is an auxiliary continuous variable and £(x) is an effective position operator of dimension N x N. This is 
defined by setting [£(x)]jj = [x]ij if (hj) ^ (N,N) and [£(x)]n,n = %n,n(x)> with 

N/2 

*^) = (jV+l)/2 ■ (A2) 

x + ie 

X + IE — . . . 

The single term xn,n{x) contains the contribution of all the states which are not occupied by the fermions. Its 
asymptotic value is Xn,n(%) — iy2N — x 2 for N — > oo. 

An expression for Tr^v G(x) is obtained from the Green's function element Qx } n(8, x) between the first and the last 
occupied state by using the expression 

— lnt/i n[o,x) = Tr?ln = TT7 In — tj 

05 85 dctlx + ^-^x^e] 85 [ x + 5 - £i(x) + ie] 

q N 1 N _ 1 

_ 95 X + S _ +i£ ~ X + 5 _ + i£ 

where ^i(x) are the eigenvalues of the operator £(x). Therefore, the particle density is given by 



nix) = lim ImTrArG(x) = — lim Im 

7T e^0+ 7T £-.0+ 
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lnt/i jA r((5,x) 



(A4) 



<5=0 



The matrix element Qi,n(5, x) can be calculated by a further renormalization |l5| of ^(x), without explicitly inverting 
the matrix (x + 5 — £(x) + ie). As in the method presented in Sec. 0, convergence is achieved by setting M = 10 7 
and £ = 0.001. 

The same procedure can be used to evaluate the momentum distribution as the trace of an N x N operator 
7r(p), defined analogously to £(x). Due to the symmetry between the operators x and p in the harmonic oscillator 
Hamiltonian, we obtain essentially the same final expression as for of the particle density nix). Thus, the momentum 
distribution for N fermions under ID harmonic confinement shows N oscillations. 

The mapping between hard-core bosons and non-interacting fermions in ID |l0| allows us to evaluate the particle 
density profile for mesoscopic systems of strongly correlated bosons by the same method. However, the mapping is 
restricted to real-space properties and does not apply to the momentum distribution. 
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(a) (b) 

FIG. 1. Exact kinetic energy density T(x) (a) and momentum flux density II(a;) (b), in units of hu/aho as functions of the 
spatial coordinate x (in units of a^o = \JhjmuS] for N = 4 fermions (dashed line), 12 fermions (dotted line) and 24 fermions 
(solid line) in a ID harmonic potential. The inset in (b) shows the derivative dH(x)/dx. 
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FIG. 3. Quantal fluctuations AT(x) of the kinetic energy density for a system of two fermions evaluated numerically by the 
Green's function method (dotted line), compared with the result of an analytic calculation (continuous line). The difference 
between the two curves is shown in the inset. The units are the same as in Fig. 0. 
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FIG. 4. Local kinetic energy density (in units of N~ 3 ^ 2 hui/aho) as a function of the particle density (in units of N^^a^l) 
for N = 12. The circles indicate the regions in which the function is multi- valued. The inset shows an enlargement of one of 
these regions, indicated by the arrow. The function T(n) is multi-valued around all local maxima of the density profile. 
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FIG. 5. Particle density n (in units of N 1 ^ 2 a ho 1 ) as a function of ni oc (in units of hoj/N) for various values of N as compared 
with the LDA profile. 
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